
########################################
############# Figure A3 ################
########################################

library(sf)
library(terra)
library(dplyr)
library(tidyr)
library(haven)
library(ggplot2)
library(exactextractr)

# nightlight rasters
nl2014 <- rast("Data/nightlights/F152014ext.tif")
nl2019 <- rast("Data/nightlights/F152019ext.tif")

# muncipalities
municips <- st_read("Data/communes_shape/communes-20181110.shp")
municips <- st_transform(municips, crs = crs(nl2019)) # ensure same crs

# zonal stats
zonal_nl2014 <- exact_extract(nl2014, municips, c('sum', 'mean'))
zonal_nl2019 <- exact_extract(nl2019, municips, c('sum', 'mean'))

# add to sf
municips$nl2014sum <- zonal_nl2014$sum
municips$nl2014mean <- zonal_nl2014$mean
municips$nl2019sum <- zonal_nl2019$sum
municips$nl2019mean <- zonal_nl2019$mean

# replace zeros
municips <- municips %>%
  mutate(
    nl2014sum = replace_na(zonal_nl2014$sum, 0),
    nl2014mean = replace_na(zonal_nl2014$mean, 0),
    nl2019sum = replace_na(zonal_nl2019$sum, 0),
    nl2019mean = replace_na(zonal_nl2019$mean, 0)
  )

## plot

continental_bbox <- st_as_sfc(st_bbox(c(xmin = -5, xmax = 10, ymin = 41, ymax = 51), crs = 4326))
# Ensure the shapefile is in WGS84 (EPSG:4326)
municips <- st_transform(municips, crs = 4326)

# Filter municipalities within the bounding box
continental_france <- municips %>%
  st_filter(continental_bbox)

# logged variables
continental_france$log_nl2019 <- log(continental_france$nl2019sum + 1) 
continental_france$log_diff_nl <- log(continental_france$nl2019sum + 1) - log(continental_france$nl2014sum + 1)

nl2019 <- ggplot(data = continental_france) +
  geom_sf(aes(fill = log_nl2019), color = NA) +
  scale_fill_viridis_c(option = "plasma", name = "Log Nightlights\n(2019)") +
  theme_minimal() +
  labs(title = "Logged Nightlight Emissions (2019)", 
       caption = "Data source: Extending the DMSP Nighttime Lights Time Series (Elvidge et al. 1997, 2021; Gosh et al. 2021") +
  theme(legend.position = "right")


# Plot Difference in Logged Nightlight Emissions (2014 to 2019)
nl2014_9 <- ggplot(data = continental_france) +
  geom_sf(aes(fill = log_diff_nl), color = NA) +
  scale_fill_viridis_c(option = "plasma", name = "Change in Log\nNightlights") +
  theme_minimal() +
  labs(title = "Change in Logged Nightlight Emissions (2014-2019)", 
       caption = "Data source: Extending the DMSP Nighttime Lights Time Series (Elvidge et al. 1997, 2021; Gosh et al. 2021") +
  theme(legend.position = "right")

# Export 
ggsave("Figures/figA3a.jpg", plot = nl2019, width = 8, height = 6)
ggsave("Figures/figA3b.jpg", plot = nl2014_9, width = 8, height = 6)



